
#helper functions---------------------------------------------------

make_simulation = function(tau, assignment_prob, alpha = alpha, interaction = NULL, int_effect = NULL){
  
  possible.ns <- seq(from=100, to=2000, by=50) # The sample sizes we'll be considering 
  powers <- rep(NA, length(possible.ns)) # Empty object to collect simulation estimates 
  alpha <- alpha # Standard significance level 
  sims <- 500 # Number of simulations to conduct for each N 
  sim = NULL
  

#### Outer loop to vary the number of subjects #### 
      for (j in 1:length(possible.ns)){ 
        
        if(interaction == "none"){
          N <- possible.ns[j] # Pick the jth value for N 
          
          Y0 <- rnorm(n=N, mean=0.43, sd=0.2) # control potential outcome 
          Y1 <- Y0 + tau # treatment potential outcome                                   
          significant.experiments <- rep(NA, sims) # Empty object to count significant experiments
          
      #### Inner loop to conduct experiments "sims" times over for each N #### 
            for (i in 1:sims){
              Z.sim <- rbinom(n=N, size=1, prob= assignment_prob) # Do a random assignment 
              Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment 
              fit.sim <- lm(Y.sim ~ Z.sim) # Do analysis (Simple regression) 
              p.value <- summary(fit.sim)$coefficients[2,4] # Extract p-values 
              significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
            }
        }
        if(interaction == "pious"){
          
          N <- possible.ns[j] # Pick the jth value for N 
                            
          significant.experiments <- rep(NA, sims) # Empty object to count significant experiments
          for (i in 1:sims){
            pious = rbinom(n=N, size=1, prob= 0.5)
            
            Y0 <- rnorm(n=N, mean=0.43, sd=0.2) # control potential outcome
            
            Y1 <- Y0 + tau + (pious*(int_effect)) # treatment potential outcome       
            
            Z.sim <- rbinom(n=N, size=1, prob= assignment_prob) # Do a random assignment 
            Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment 
            fit.sim <- lm(Y.sim ~ Z.sim*pious) # Do analysis (Simple regression) 
            p.value <- try(summary(fit.sim)$coefficients[4,4], TRUE) # Extract p-values 
            if(class(p.value) == "try-error"){next}
            significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
          }
        }
        
        if(interaction == "islam"){
          
          N <- possible.ns[j] # Pick the jth value for N 
                      
          significant.experiments <- rep(NA, sims) # Empty object to count significant experiments
          for (i in 1:sims){
            islam = rbinom(n=N, size=1, prob= 0.88)
            
            Y0 <- rnorm(n=N, mean=0.43, sd=0.2) # control potential outcome
            
            Y1 <- Y0 + tau + (islam*(int_effect)) # treatment potential outcome       
            
            Z.sim <- rbinom(n=N, size=1, prob= assignment_prob) # Do a random assignment 
            Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment 
            fit.sim <- lm(Y.sim ~ Z.sim*islam) # Do analysis (Simple regression) 
            p.value <- try(summary(fit.sim)$coefficients[4,4], TRUE) # Extract p-values 
            if(class(p.value) == "try-error"){next}
            significant.experiments[i] <- (p.value <= alpha)
          }
        }
        
        if(interaction == "noncompliance"){
          
          N <- possible.ns[j] # Pick the jth value for N 
          
          significant.experiments <- rep(NA, sims) # Empty object to count significant experiments
          for (i in 1:sims){
            noncompliance = rbinom(n=N, size=1, prob= 0.42)
            
            Y0 <- rnorm(n=N, mean=0.43, sd=0.2) # control potential outcome
            
            Y1 <- Y0 + tau + (noncompliance*(int_effect)) # treatment potential outcome       
            
            Z.sim <- rbinom(n=N, size=1, prob= assignment_prob) # Do a random assignment 
            Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment 
            fit.sim <- lm(Y.sim ~ Z.sim*noncompliance) # Do analysis (Simple regression) 
            p.value <- try(summary(fit.sim)$coefficients[4,4], TRUE) # Extract p-values 
            if(class(p.value) == "try-error"){next}
            significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
          }
        }
        
        if(interaction == "jokowi"){
          
          N <- possible.ns[j] # Pick the jth value for N 
          
          significant.experiments <- rep(NA, sims) # Empty object to count significant experiments
          for (i in 1:sims){
            jokowi = rbinom(n=N, size=1, prob= 0.55)
            
            Y0 <- rnorm(n=N, mean=0.43, sd=0.2) # control potential outcome
            
            Y1 <- Y0 + tau + (jokowi*(int_effect)) # treatment potential outcome       
            
            Z.sim <- rbinom(n=N, size=1, prob= assignment_prob) # Do a random assignment 
            Y.sim <- Y1*Z.sim + Y0*(1-Z.sim) # Reveal outcomes according to assignment 
            fit.sim <- lm(Y.sim ~ Z.sim*jokowi) # Do analysis (Simple regression) 
            p.value <- try(summary(fit.sim)$coefficients[4,4], TRUE) # Extract p-values 
            if(class(p.value) == "try-error"){next}
            significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
          }
        }
        
        if(interaction == "none"){
        sim = 
          sim %>%
          bind_rows(c(out = mean(significant.experiments, na.rm = T), N = possible.ns[j], tau = tau))
        } else if(interaction %in% c("islam", "pious", "noncompliance", "jokowi")){
          sim = 
            sim %>%
            bind_rows(c(out = mean(significant.experiments, na.rm = T), N = possible.ns[j], tau = int_effect))
          
        }
        
        
      }
  
  return(sim)
}

#second helper function for making the graphs

make_plot = function(df, num_obs, interaction = NULL){
  
  if(interaction == "none"){lab = expression(beta[1])}
  if(interaction %in% c("pious", "islam", "noncompliance", "jokowi")){lab = expression(beta[3])}
  
  df %>%
  ggplot(aes(x = N, y = out*100, color = factor(tau))) +
    geom_point(size = 1) +
#    geom_line() +
    geom_vline(xintercept = num_obs, color = "red", linetype = "dashed") +
    geom_hline(yintercept = 95, color = "red", linetype = "dashed") +
    theme_bw() +
    ylab("Probability of Detecting Effect (%)") +
    xlab("Number of Observations") +
    labs(color=lab) +
    scale_color_viridis_d() +
    theme(text=element_text(size=12, family="Times")) + 
    theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = .2)) + 
    theme(strip.background =element_rect(fill="white")) +
    theme(panel.border = element_blank(),
          axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
          axis.line.y = element_line(size = 0.5, linetype = "solid", colour = "black")) %>%
    return()
}


#H1 power calculation----------------------------------------------------
H1 = 
  make_simulation(tau = 0.05, assignment_prob = 0.8, alpha = 0.05, interaction = "none") %>%
  bind_rows(., make_simulation(tau = 0.06, assignment_prob = 0.8, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.08, assignment_prob = 0.8, alpha = 0.05, interaction = "none"))


H1_plot = make_plot(H1, num_obs = 1600, interaction = "none")

#H2 power calculation----------------------------------------------------

H2 = 
  make_simulation(tau = 0.05, assignment_prob = 0.5, alpha = 0.05, interaction = "none") %>%
  bind_rows(., make_simulation(tau = 0.06, assignment_prob = 0.5, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.07, assignment_prob = 0.5, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.08, assignment_prob = 0.5, alpha = 0.05, interaction = "none"))


H2_plot = make_plot(H2, num_obs = 1200, interaction = "none")

#H3 power calculation----------------------------------------------------

H3 = 
  make_simulation(tau = 0.05, assignment_prob = 0.5, alpha = 0.05, interaction = "none") %>%
  bind_rows(., make_simulation(tau = 0.06, assignment_prob = 0.5, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.07, assignment_prob = 0.5, alpha = 0.05, interaction = "none")) %>%
  bind_rows(., make_simulation(tau = 0.08, assignment_prob = 0.5, alpha = 0.05, interaction = "none"))


H3_plot = make_plot(H3, num_obs = 1200, interaction = "none")

#H4 power calculation----------------------------------------------------

H4_05 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.05)
H4_06 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.06)
H4_07 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.07)
H4_08 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.08)
H4_09 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.09)
H4_10 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "pious", int_effect = 0.10)

H4 =
  bind_rows(H4_05, H4_06) %>%
  bind_rows(., H4_07) %>%
  bind_rows(., H4_08) %>%
  bind_rows(., H4_09) %>%
  bind_rows(., H4_10)


H4_plot = make_plot(H4, num_obs = 1600, interaction = "pious")


#H5 power calculation----------------------------------------------------


H5_10 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.10)
H5_11 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.11)
H5_12 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.12)
H5_13 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.13)
H5_14 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.14)
H5_15 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "islam", int_effect = 0.15)

H5 =
  bind_rows(H5_10, H5_11) %>%
  bind_rows(., H5_12) %>%
  bind_rows(., H5_13) %>%
  bind_rows(., H5_14) %>%
  bind_rows(., H5_15)


H5_plot = make_plot(H5, num_obs = 1600, interaction = "islam")


#H6 power calculation----------------------------------------------------


H6_12 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "noncompliance", int_effect = 0.07)
H6_13 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "noncompliance", int_effect = 0.08)
H6_14 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "noncompliance", int_effect = 0.09)
H6_15 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "noncompliance", int_effect = 0.10)

H6 =
  bind_rows(H6_12, H6_13) %>%
  bind_rows(., H6_14) %>%
  bind_rows(., H6_15)


H6_plot = make_plot(H6, num_obs = 1600, interaction = "noncompliance")


#H7 power calculation----------------------------------------------------


H7_07 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "jokowi", int_effect = 0.07)
H7_08 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "jokowi", int_effect = 0.08)
H7_09 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "jokowi", int_effect = 0.09)
H7_10 = make_simulation(tau = 0.07, assignment_prob = 0.8, alpha = 0.05, interaction = "jokowi", int_effect = 0.10)

H7 =
  bind_rows(H7_07, H7_08) %>%
  bind_rows(., H7_09) %>%
  bind_rows(., H7_10)


H7_plot = make_plot(H7, num_obs = 1600, interaction = "jokowi")



#saving plots------------------------------------------------------------

ggsave(plot = H1_plot, "./outputs/sa_figure_1a.pdf", width = 5, height = 3)
ggsave(plot = H2_plot, "./outputs/sa_figure_1b.pdf", width = 5, height = 3)
ggsave(plot = H3_plot, "./outputs/sa_figure_1c.pdf", width = 5, height = 3)
ggsave(plot = H4_plot, "./outputs/sa_figure_1d.pdf", width = 5, height = 3)
ggsave(plot = H5_plot, "./outputs/sa_figure_1e.pdf", width = 5, height = 3)
ggsave(plot = H6_plot, "./outputs/sa_figure_1f.pdf", width = 5, height = 3)
ggsave(plot = H7_plot, "./outputs/sa_figure_1g.pdf", width = 5, height = 3)



